{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "dl_implement_av.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.8"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3pSr3puhY6oC"
      },
      "source": [
        "Reading and pre-processing of time-series data for 4 stations"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZyT4BH9Db7FX"
      },
      "source": [
        "                                                                                                                                                                                                                                                                        from math import sqrt\n",
        "import numpy as np\n",
        "from numpy import concatenate\n",
        "from pandas import read_csv\n",
        "from pandas import DataFrame\n",
        "from pandas import concat\n",
        "from sklearn.preprocessing import MinMaxScaler\n",
        "from sklearn.preprocessing import LabelEncoder\n",
        "from sklearn.metrics import mean_squared_error\n",
        "from keras.models import Sequential\n",
        "from keras.layers import Dense\n",
        "from keras.layers import LSTM\n",
        "from pandas import read_csv\n",
        "from datetime import datetime\n",
        "from keras.layers import Bidirectional\n",
        "from keras.layers import Dropout\n",
        "import datetime\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "import time\n",
        "import os\n",
        "import seaborn as sn\n",
        "import copy\n",
        "## Animesh's Update\n",
        "\n",
        "# load data\n",
        "def parse(x):\n",
        "\treturn datetime.strptime(x, '%Y %m %d %H')\n",
        "\n",
        "# Anand Vihar \n",
        "dfav = read_csv('data_new/anand_vihar_jan19_dec20_8h.csv')\n",
        "dfav.columns = ['From Date','PM10','NO','NO2','WS','CO','Benzene','NOx','Ozone','SO2','NH3','Toluene','PM2.5']\n",
        "\n",
        "# Bawana\n",
        "dfbw= read_csv('data_new/Bawana_jan19-dec20_8h.csv')\n",
        "dfbw.columns = ['From Date','PM10','NO','NO2','WS','CO','Benzene','NOx','Ozone','SO2','NH3','Toluene','PM2.5']\n",
        "\n",
        "# Vivek Vihar\n",
        "dfvv = read_csv('data_new/vivek_vihar_jan19_dec20_8h.csv')\n",
        "dfvv.columns = ['From Date','PM10','NO','NO2','WS','CO','Benzene','NOx','Ozone','SO2','NH3','Toluene','PM2.5']\n",
        "\n",
        "\n",
        "dfwz = read_csv('data_new/Wazirpur_8h_jan19-dec20.csv')\n",
        "dfwz.columns = ['From Date','PM10','NO','NO2','WS','CO','Benzene','NOx','Ozone','SO2','NH3','Toluene','PM2.5']\n",
        "\n",
        "\n",
        "# DTU\n",
        "dfdt = read_csv('data_new/DTU_8h_jan19-dec20.csv')\n",
        "dfdt.columns = ['From Date','PM10','NO','NO2','WS','CO','Benzene','NOx','Ozone','SO2','NH3','Toluene','PM2.5']\n",
        "\n",
        "\n",
        "\n",
        "# Edit it accordingly for 4 different stations: Currently set-up for Anand Vihar\n",
        "df= dfav\n",
        "\n",
        "print(\"Shape of df:\",df.shape)\n",
        "\n",
        "\n",
        "print(df.head(5))\n",
        "\n",
        "\n",
        "df = df.replace(to_replace='None', value=np.nan)\n",
        "df = df.dropna(axis = 0, how = 'all')\n",
        "df['PM10'].fillna(df['PM10'].median(), inplace=True)\n",
        "df['NO'].fillna(df['NO'].median(), inplace=True)\n",
        "df['NO2'].fillna(df['NO2'].median(), inplace=True)\n",
        "df['WS'].fillna(df['WS'].median(), inplace=True)\n",
        "df['CO'].fillna(df['CO'].median(), inplace=True)\n",
        "df['Benzene'].fillna(df['Benzene'].median(), inplace=True)\n",
        "df['NOx'].fillna(df['NOx'].median(), inplace=True)\n",
        "df['Ozone'].fillna(df['Ozone'].median(), inplace=True)\n",
        "df['SO2'].fillna(df['SO2'].median(), inplace=True)\n",
        "df['NH3'].fillna(df['NH3'].median(), inplace=True)\n",
        "df['Toluene'].fillna(df['Toluene'].median(), inplace=True)\n",
        "df['PM2.5'].fillna(df['PM2.5'].median(), inplace=True)\n",
        "\n",
        "print(\" \\nCount total NaN at each column in a DataFrame : \\n\\n\", df.isnull().sum())\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EsQ5K-eUyrZJ"
      },
      "source": [
        "df.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b30HPz1R4TWN"
      },
      "source": [
        "features = df.loc[:, df.columns != 'From Date']\r\n",
        "features = features.astype(float)\r\n",
        "features.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GkFhZ10wsP97"
      },
      "source": [
        "# Table for range, mean and std for different features\r\n",
        "feat = features['PM10'].values\r\n",
        "print(\"PM10 :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['NO'].values\r\n",
        "print(\"NO :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['NO2'].values\r\n",
        "print(\"NO2 :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['WS'].values\r\n",
        "print(\"WS :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['CO'].values\r\n",
        "print(\"CO :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['Benzene'].values\r\n",
        "print(\"Benzene :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['NOx'].values\r\n",
        "print(\"NOx :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['Ozone'].values\r\n",
        "print(\"Ozone :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['SO2'].values\r\n",
        "print(\"SO2 :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['NH3'].values\r\n",
        "print(\"NH3 :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['Toluene'].values\r\n",
        "print(\"Toluene :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))\r\n",
        "feat = features['PM2.5'].values\r\n",
        "print(\"PM2.5 :\", \"Min = \",min(feat), \", Max = \", max(feat) , \", Mean = \", np.mean(feat), \", Std = \",np.std(feat))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N9wYS3lsZNmi"
      },
      "source": [
        "Correlation Matrix and Heatmap for features"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2fS1J3BZ_cQs"
      },
      "source": [
        "\r\n",
        "corr_matrix = features.corr()\r\n",
        "f, ax = plt.subplots(figsize=(11, 15))\r\n",
        "heatmap = sn.heatmap(corr_matrix,\r\n",
        "                      #mask = mask,\r\n",
        "                      square = True,\r\n",
        "                      linewidths = .5,\r\n",
        "                      cmap = 'coolwarm',\r\n",
        "                      cbar_kws = {'shrink': .4,\r\n",
        "                                'ticks' : [-1, -.5, 0, 0.5, 1]},\r\n",
        "                      vmin = -1,\r\n",
        "                      vmax = 1,\r\n",
        "                      annot = True,\r\n",
        "                      annot_kws = {\"size\": 12})\r\n",
        "\r\n",
        "#add the column names as labels\r\n",
        "ax.set_yticklabels(corr_matrix.columns, rotation = 0)\r\n",
        "ax.set_xticklabels(corr_matrix.columns)\r\n",
        "\r\n",
        "sn.set_style({'xtick.bottom': True}, {'ytick.left': True})"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qcp7Vu_kZhL6"
      },
      "source": [
        "Preparation of dataframe for analysis"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8cNGumhQb7Fi"
      },
      "source": [
        "df['Date'] = pd.to_datetime(df['From Date']) \n",
        "df['year'] = df['Date'].dt.year \n",
        "df['month'] = df['Date'].dt.month \n",
        "df['day'] = df['Date'].dt.day \n",
        "df['hour'] = df['Date'].dt.hour \n",
        "df['minute'] = df['Date'].dt.minute\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FPWoN7Xenms0"
      },
      "source": [
        "df.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1k157kZ2ZmjO"
      },
      "source": [
        "Splitting into train/test for different LSTM models depending on months covering pre-lockdown and during/post-lockdown respectively."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Rl6pRuIEb7Fj"
      },
      "source": [
        "df_train = df[(df[\"year\"] <2020) | (df['month']<=5)]\n",
        "\n",
        "df_test=df[(df[\"year\"] >2019) & (df['month']<=12) & (df['month']>5)]\n",
        "df.to_csv('pollution.csv')\n",
        "df_train.to_csv('train_pollution.csv')\n",
        "df_test.to_csv('test_pollution.csv')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Rb6kKBbsZ-F5"
      },
      "source": [
        "For season-wise training and testing. Uncomment this for corresponding seasonal analysis."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Tz2NZR9yTB3a"
      },
      "source": [
        "## For season-wise training and testing\r\n",
        "\"\"\"\r\n",
        "df_train = df[((df[\"year\"] == 2019) ) & (df['month'] >= 2) & (df['month'] <=9)]\r\n",
        "df_test=df[(df[\"year\"]==2020) & (df['month'] >= 2) & (df['month'] <=9)]\r\n",
        "df.to_csv('pollution.csv')\r\n",
        "df_train.to_csv('train_pollution.csv')\r\n",
        "df_test.to_csv('test_pollution.csv')\"\"\""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PHlXsMmHe8Yt"
      },
      "source": [
        "# Note: df_july has to be used only for seasonal analysis\r\n",
        "df_july = df[((df[\"year\"] == 2020) ) & (df['month'] == 9)  ]\r\n",
        "#df_july = df[((df[\"year\"] == 2020) ) & (df['month'] >= 8) & (df['month'] <=9)]\r\n",
        "df_july.to_csv('july_test.csv')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wS1NvpsWbGEv"
      },
      "source": [
        "Plotting the time-series data for train"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rj_Pw1blD04c"
      },
      "source": [
        "x_val = np.linspace(1,df_train.shape[0],num = df_train.shape[0])\r\n",
        "plot_list = np.array(df_train['PM2.5'].values)\r\n",
        "y_val = np.array([(float)(k) for k in plot_list])\r\n",
        "fig = plt.figure(figsize = (12,12))\r\n",
        "ax = fig.add_subplot(111) \r\n",
        "plt.tick_params(labelsize=25)\r\n",
        "params = {'legend.fontsize': 25, 'legend.handlelength': 2}\r\n",
        "plt.rcParams.update(params)\r\n",
        "surrogate_plot = ax.plot(x_val,y_val,linestyle='-', linewidth= 1, color= 'b', label= 'PM2.5')\r\n",
        "ax.set_xlabel('Time (Samples)',size= 25)\r\n",
        "ax.set_ylabel('PM2.5 Value', size= 25)\r\n",
        "ax.set_xlim([0,np.amax(x_val)])   \r\n",
        "ax.legend(loc='best')\r\n",
        "fig.tight_layout()\r\n",
        "fig.subplots_adjust(top=0.88)\r\n",
        "plt.show()\r\n",
        "plt.clf()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t5-1s5egbKcC"
      },
      "source": [
        "Plotting the time-series data for test"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YM3HmkdhD0xs"
      },
      "source": [
        "x_val = np.linspace(1,df_test.shape[0],num = df_test.shape[0])\r\n",
        "plot_list = np.array(df_test['PM2.5'].values)\r\n",
        "y_val = np.array([(float)(k) for k in plot_list])\r\n",
        "fig = plt.figure(figsize = (12,12))\r\n",
        "ax = fig.add_subplot(111) \r\n",
        "plt.tick_params(labelsize=25)\r\n",
        "params = {'legend.fontsize': 25, 'legend.handlelength': 2}\r\n",
        "plt.rcParams.update(params)\r\n",
        "surrogate_plot = ax.plot(x_val,y_val,linestyle='-', linewidth= 1, color= 'r', label= 'PM2.5')\r\n",
        "ax.set_xlabel('Time Steps',size= 25)\r\n",
        "ax.set_ylabel('PM2.5 Concentration', size= 25)\r\n",
        "ax.set_xlim([0,np.amax(x_val)])   \r\n",
        "ax.legend(loc='best')\r\n",
        "fig.tight_layout()\r\n",
        "fig.subplots_adjust(top=0.88)\r\n",
        "plt.show()\r\n",
        "plt.clf()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Mw8ADvNYb7Fk"
      },
      "source": [
        "dataset_train = read_csv('train_pollution.csv', header=0, index_col=0)\n",
        "dataset_test = read_csv('test_pollution.csv', header=0, index_col=0)\n",
        "values_train = dataset_train.values\n",
        "values_test = dataset_test.values\n",
        "values_train=values_train[:,1:-6]\n",
        "values_test=values_test[:,1:-6]\n",
        "print(values_train.shape,values_test.shape)\n",
        "n_hours=1\n",
        "n_out= 1\n",
        "dataset = read_csv('pollution.csv', header=0, index_col=0)\n",
        "values = dataset.values\n",
        "print(values.shape)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kER_lziobR09"
      },
      "source": [
        "Defining different models and trial experiments"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4F6R8UlBb7Fo"
      },
      "source": [
        "from keras.layers import RepeatVector\n",
        "from keras.layers import TimeDistributed\n",
        "from numpy import array\n",
        "import datetime\n",
        "\n",
        "def split_sequences(sequences, n_steps_in, n_steps_out):\n",
        "    X, y = list(), list()\n",
        "    for i in range(len(sequences)):\n",
        "    # find the end of this pattern\n",
        "        end_ix = i + n_steps_in\n",
        "        out_end_ix = end_ix + n_steps_out-1\n",
        "    # check if we are beyond the dataset\n",
        "        if out_end_ix > len(sequences):\n",
        "            break\n",
        "        # gather input and output parts of the pattern\n",
        "        \n",
        "        seq_x, seq_y = sequences[i:end_ix, :-1], sequences[end_ix-1:out_end_ix, -1]\n",
        "        X.append(seq_x)\n",
        "        y.append(seq_y)\n",
        "    return array(X), array(y)\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "def MODEL_FNN_adam(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv):\n",
        "    print(train_X.shape)\n",
        "    train_X=np.reshape(train_X,(train_X.shape[0],n_hours*n_features))\n",
        "    test_X=np.reshape(test_X,(test_X.shape[0],n_hours*n_features))\n",
        "    val_X=np.reshape(val_X,(val_X.shape[0],n_hours*n_features))\n",
        "\n",
        "    model = Sequential()\n",
        "    model.add(Dense(64, input_dim=train_X.shape[1], activation='relu'))\n",
        "    model.add(Dense(32, activation='relu'))\n",
        "    model.add(Dense(n_out))\n",
        "    model.compile(loss='mse', optimizer='adam')\n",
        "    # fit the keras model on the dataset\n",
        "    model.fit(train_X, train_y, epochs=200, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle=False)\n",
        "    yhat11 = model.predict(test_X, verbose=0)\n",
        "    yhat11_inv = ((maxv - minv) * yhat11  + minv)\n",
        "    test_y_inv = ((maxv - minv) * test_y  + minv)\n",
        "    print(\"Predicted = \", yhat11_inv)\n",
        "    print(\"Actual = \", test_y_inv)\n",
        "    err = np.sqrt(((np.sum((yhat11_inv-test_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step is =\",err )\n",
        "    yhat12 = model.predict(train_X, verbose = 0)\n",
        "    yhat12_inv = ((maxv - minv) * yhat12  + minv)\n",
        "    train_y_inv = ((maxv - minv) * train_y  + minv)\n",
        "    err2 = np.sqrt(((np.sum((yhat12_inv-train_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step for train is =\",err2 )\n",
        "    return err , err2\n",
        "\n",
        "\n",
        "\n",
        "def MODEL_LSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv):\n",
        "    model = Sequential()\n",
        "    model.add(LSTM(50, activation='relu', input_shape=(n_hours,n_features)))\n",
        "    model.add(Dense(n_out))\n",
        "    model.compile(optimizer='adam', loss='mse')\n",
        "# fit model\n",
        "    model.fit(train_X, train_y, epochs=200, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle=False)\n",
        "    yhat11 = model.predict(test_X, verbose=0)\n",
        "    yhat11_inv = ((maxv - minv) * yhat11  + minv)\n",
        "    test_y_inv = ((maxv - minv) * test_y  + minv)\n",
        "    print(\"Predicted = \", yhat11_inv)\n",
        "    print(\"Actual = \", test_y_inv)\n",
        "    err = np.sqrt(((np.sum((yhat11_inv-test_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step is =\",err )\n",
        "    yhat12 = model.predict(train_X, verbose = 0)\n",
        "    yhat12_inv = ((maxv - minv) * yhat12  + minv)\n",
        "    train_y_inv = ((maxv - minv) * train_y  + minv)\n",
        "    err2 = np.sqrt(((np.sum((yhat12_inv-train_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step for train is =\",err2 )\n",
        "    return err , err2\n",
        "\n",
        "def MODEL_Bi_LSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv):\n",
        "    model = Sequential()\n",
        "    model.add(Bidirectional(LSTM(50, activation='relu',input_shape=(n_hours,n_features))))\n",
        "    model.add(Dense(n_out))\n",
        "    model.compile(optimizer='adam', loss='mse')\n",
        "# fit model\n",
        "    # Important note : epochs = 50 for seasonal model\n",
        "    model.fit(train_X, train_y, epochs=200, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle=False)\n",
        "    yhat11 = model.predict(test_X, verbose=0)\n",
        "    yhat11_inv = ((maxv - minv) * yhat11  + minv)\n",
        "    test_y_inv = ((maxv - minv) * test_y  + minv)\n",
        "    print(\"Predicted = \", yhat11_inv)\n",
        "    print(\"Actual = \", test_y_inv)\n",
        "    err = np.sqrt(((np.sum((yhat11_inv-test_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step is =\",err )\n",
        "    yhat12 = model.predict(train_X, verbose = 0)\n",
        "    yhat12_inv = ((maxv - minv) * yhat12  + minv)\n",
        "    train_y_inv = ((maxv - minv) * train_y  + minv)\n",
        "    err2 = np.sqrt(((np.sum((yhat12_inv-train_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step for train is =\",err2 )\n",
        "    return err , err2\n",
        "\n",
        "\n",
        "def MODEL_EN_DC(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv):\n",
        "    train_y = train_y.reshape((train_y.shape[0], train_y.shape[1], 1))\n",
        "    test_y = test_y.reshape((test_y.shape[0], test_y.shape[1], 1))\n",
        "    val_y = val_y.reshape((val_y.shape[0], val_y.shape[1], 1))\n",
        "    model = Sequential()\n",
        "    model.add(LSTM(50, activation='relu',input_shape=(n_hours,n_features)))\n",
        "    model.add(RepeatVector(n_out))\n",
        "    model.add(LSTM(50, activation='relu', return_sequences=True))\n",
        "    model.add(TimeDistributed(Dense(1)))\n",
        "    model.compile(optimizer='adam', loss='mse')\n",
        "# fit model\n",
        "    model.fit(train_X, train_y, epochs=200, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle=False)\n",
        "    yhat11 = model.predict(test_X, verbose=0)\n",
        "    yhat11_inv = ((maxv - minv) * yhat11  + minv)\n",
        "    test_y_inv = ((maxv - minv) * test_y  + minv)\n",
        "    print(\"Predicted = \", yhat11_inv)\n",
        "    print(\"Actual = \", test_y_inv)\n",
        "    err = np.sqrt(((np.sum((yhat11_inv-test_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step is =\",err )\n",
        "    yhat12 = model.predict(train_X, verbose = 0)\n",
        "    yhat12_inv = ((maxv - minv) * yhat12  + minv)\n",
        "    train_y_inv = ((maxv - minv) * train_y  + minv)\n",
        "    err2 = np.sqrt(((np.sum((yhat12_inv-train_y_inv)**2,axis=0))))\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step for train is =\",err2 )\n",
        "    return err , err2\n",
        "\n",
        "           \n",
        "            \n",
        "            \n",
        "            \n",
        "n_out_max,n_hour_max= 10,5\n",
        "out_range=[10]\n",
        "dataset_train = read_csv('train_pollution.csv', header=0, index_col=0)\n",
        "dataset_test = read_csv('test_pollution.csv', header=0, index_col=0)\n",
        "values_train = dataset_train.values\n",
        "values_test = dataset_test.values\n",
        "values_train=values_train[:,1:-6]\n",
        "values_test=values_test[:,1:-6]\n",
        "\n",
        "values_train = values_train.astype('float32')\n",
        "values_test = values_test.astype('float32')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QUY8BnENcKur"
      },
      "source": [
        "Data preparation for training/testing of models"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HnCO4I-bGmhH"
      },
      "source": [
        "# Run this cell for general analysis of models\r\n",
        "dataset_train = read_csv('train_pollution.csv', header=0, index_col=0)\r\n",
        "dataset_test = read_csv('test_pollution.csv', header=0, index_col=0)\r\n",
        "values_train = dataset_train.values\r\n",
        "values_test = dataset_test.values\r\n",
        "values_train=values_train[:,1:-6]\r\n",
        "values_test=values_test[:,1:-6]\r\n",
        "\r\n",
        "values_train = values_train.astype('float32')\r\n",
        "values_test = values_test.astype('float32')\r\n",
        "scaler = MinMaxScaler(feature_range=(0, 1))\r\n",
        "concat = np.concatenate((values_train,values_test),axis = 0)\r\n",
        "maxv = max(concat[:,-1])\r\n",
        "minv = min(concat[:,-1])\r\n",
        "scaled = scaler.fit_transform(concat)\r\n",
        "scaled_train = scaled[:values_train.shape[0],:]\r\n",
        "scaled_test = scaled[values_train.shape[0]:,:]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NjeCmzukcFJd"
      },
      "source": [
        "Data preparation here is only for 1 month comparison for September using SBD-LSTM\r\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "v-L19iC6fvIU"
      },
      "source": [
        "# Run this cell only for seasonal analysis of September\r\n",
        "data_july = read_csv('july_test.csv', header=0, index_col=0)\r\n",
        "values_july = data_july.values\r\n",
        "values_july = values_july[:,1:-6]\r\n",
        "values_july = values_july.astype('float32')\r\n",
        "scaler = MinMaxScaler(feature_range=(0, 1))\r\n",
        "concat = np.concatenate((values_train,values_test,values_july),axis = 0)\r\n",
        "maxv = max(concat[:,-1])\r\n",
        "minv = min(concat[:,-1])\r\n",
        "scaled = scaler.fit_transform(concat)\r\n",
        "scaled_train = scaled[:values_train.shape[0],:]\r\n",
        "scaled_test = scaled[values_train.shape[0]: (values_train.shape[0] + values_test.shape[0]),:]\r\n",
        "scaled_july = scaled[(values_train.shape[0] + values_test.shape[0]):,:]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-_pSdwG7df3P"
      },
      "source": [
        "Training the FNN model using Adam"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ppNh91-Cb7F4"
      },
      "source": [
        "time1 = time.time()\n",
        "outfile_adam = open('results/AnandVihar/results_adam.txt','a+')\n",
        "train_adam = open('results/AnandVihar/adam_train.txt','a+')\n",
        "for i,k in enumerate(out_range):\n",
        "    for j in range(1,31):\n",
        "        for m in range(5,6):\n",
        "            n_out=k\n",
        "            n_hours=m\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\n",
        "            \n",
        "            train_X, train_y = split_sequences(scaled_train,n_hours,n_out)\n",
        "            val_X, val_y = split_sequences(scaled_test,n_hours,n_out)\n",
        "            print(val_X.shape,val_y.shape)\n",
        "            test_X, test_y = val_X[288:,:,:],val_y[288:,:]\n",
        "            val_X, val_y = val_X[0:288,:,:],val_y[0:288,:]\n",
        "            print(\"train_X.shape, train_y.shape, test_X.shape, test_y.shape=\",train_X.shape, train_y.shape, test_X.shape, test_y.shape)\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\n",
        "            n_features = train_X.shape[2]\n",
        "            err , err2 = MODEL_FNN_adam(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv)\n",
        "            allres = np.array(err)\n",
        "            allres2 = np.array(err2)\n",
        "            for item in allres:\n",
        "              item = round(item,4)\n",
        "              outfile_adam.write(str(item) + \"  \")\n",
        "            for item in allres2:\n",
        "              item = round(item,4)\n",
        "              train_adam.write(str(item) + \"  \")  \n",
        "            #np.savetxt(outfile_adam,  allres  , fmt='%1.4f', newline=' ' )\n",
        "            #np.savetxt(train_adam,  allres2  , fmt='%1.4f', newline=' ' )\n",
        "        outfile_adam.write(\"\\n\")    \n",
        "        train_adam.write(\"\\n\")   \n",
        "\n",
        "outfile_adam.close()\n",
        "train_adam.close()\n",
        "time2 = time.time()        \n",
        "print(\"MODEL_FNN_adam:\",time2-time1)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KcW4NkyPdn23"
      },
      "source": [
        "Training LSTM model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JYQ1BnZib7F5"
      },
      "source": [
        "time1 = time.time()\n",
        "outfile_lstm = open('results/AnandVihar/results_lstm.txt','a+')\n",
        "train_lstm = open('results/AnandVihar/lstm_train.txt','a+')\n",
        "for i,k in enumerate(out_range):\n",
        "    for j in range(1,31):\n",
        "        for m in range(5,6):\n",
        "            n_out=k\n",
        "            n_hours=m\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\n",
        "            \n",
        "            train_X, train_y = split_sequences(scaled_train,n_hours,n_out)\n",
        "            val_X, val_y = split_sequences(scaled_test,n_hours,n_out)\n",
        "            print(val_X.shape,val_y.shape)\n",
        "            test_X, test_y = val_X[288:,:,:],val_y[288:,:]\n",
        "            val_X, val_y = val_X[0:288,:,:],val_y[0:288,:]\n",
        "            print(\"train_X.shape, train_y.shape, test_X.shape, test_y.shape=\",train_X.shape, train_y.shape, test_X.shape, test_y.shape)\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\n",
        "            n_features = train_X.shape[2]\n",
        "            err , err2=MODEL_LSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv)\n",
        "            allres = np.array(err)\n",
        "            allres2 = np.array(err2)\n",
        "            for item in allres:\n",
        "              item = round(item,4)\n",
        "              outfile_lstm.write(str(item) + \"  \")\n",
        "            for item in allres2:\n",
        "              item = round(item,4)\n",
        "              train_lstm.write(str(item) + \"  \")\n",
        "            #np.savetxt(outfile_lstm,  allres  , fmt='%1.4f', newline=' ' )\n",
        "            #np.savetxt(train_lstm,  allres2  , fmt='%1.4f', newline=' ' )\n",
        "        outfile_lstm.write(\"\\n\")    \n",
        "        train_lstm.write(\"\\n\") \n",
        "\n",
        "outfile_lstm.close()    \n",
        "train_lstm.close()             \n",
        "time2 = time.time()        \n",
        "print(\"MODEL_LSTM:\",time2-time1)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EPWzltt1duF0"
      },
      "source": [
        "Training Bi-LSTM model\r\n",
        "\r\n",
        "> Note: For seasonal model the split interval for validation split will be after 365 data points.\r\n",
        "\r\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jVSXg1SQb7F5"
      },
      "source": [
        "outfile_bi_lstm = open('results/AnandVihar/results_bi_lstm.txt','a+')\n",
        "train_bi_lstm = open('results/AnandVihar/bi_lstm_train.txt','a+')\n",
        "time1 = time.time()\n",
        "for i,k in enumerate(out_range):\n",
        "    for j in range(1,31):\n",
        "        for m in range(5,6):\n",
        "            n_out=k\n",
        "            n_hours=m\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\n",
        "            \n",
        "            train_X, train_y = split_sequences(scaled_train,n_hours,n_out)\n",
        "            val_X, val_y = split_sequences(scaled_test,n_hours,n_out)\n",
        "            print(val_X.shape,val_y.shape)\n",
        "            test_X, test_y = val_X[288:,:,:],val_y[288:,:]\n",
        "            val_X, val_y = val_X[0:288,:,:],val_y[0:288,:]\n",
        "            \"\"\"\n",
        "            Uncomment below two lines for Seasonal Analysis and comment above two lines\n",
        "            \"\"\"\n",
        "            #test_X, test_y = val_X[365:,:,:],val_y[365:,:]\n",
        "            #val_X, val_y = val_X[0:365,:,:],val_y[0:365,:]\n",
        "            print(\"train_X.shape, train_y.shape, test_X.shape, test_y.shape=\",train_X.shape, train_y.shape, test_X.shape, test_y.shape)\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\n",
        "            n_features = train_X.shape[2]\n",
        "            err , err2=MODEL_Bi_LSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv)\n",
        "            allres = np.array(err)\n",
        "            allres2 = np.array(err2)\n",
        "            for item in allres:\n",
        "              item = round(item,4)\n",
        "              outfile_bi_lstm.write(str(item) + \"  \")\n",
        "            for item in allres2:\n",
        "              item = round(item,4)\n",
        "              train_bi_lstm.write(str(item) + \"  \")\n",
        "            #np.savetxt(outfile_bi_lstm,  allres  , fmt='%1.4f', newline=' ' )\n",
        "            #np.savetxt(train_bi_lstm,  allres2  , fmt='%1.4f', newline=' ' )\n",
        "        outfile_bi_lstm.write(\"\\n\")    \n",
        "        train_bi_lstm.write(\"\\n\") \n",
        "            \n",
        "outfile_bi_lstm.close()    \n",
        "train_bi_lstm.close()    \n",
        "time2 = time.time()\n",
        "print(\"MODEL_Bi_LSTM:\",time2-time1)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A9p6cqFXemSh"
      },
      "source": [
        "Training Encoder-Decoder LSTM model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bQIYp4UCb7F6"
      },
      "source": [
        "outfile_en_dc = open('results/AnandVihar/results_en_dc.txt','a+')\n",
        "train_en_dc = open('results/AnandVihar/en_dc_train.txt','a+')\n",
        "time1 = time.time()\n",
        "out_range=[10]\n",
        "for i,k in enumerate(out_range):\n",
        "    for j in range(1,31):\n",
        "        for m in range(5,6):\n",
        "            n_out=k\n",
        "            n_hours=m\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\n",
        "            \n",
        "            train_X, train_y = split_sequences(scaled_train,n_hours,n_out)\n",
        "            val_X, val_y = split_sequences(scaled_test,n_hours,n_out)\n",
        "            print(val_X.shape,val_y.shape)\n",
        "            test_X, test_y = val_X[288:,:,:],val_y[288:,:]\n",
        "            val_X, val_y = val_X[0:288,:,:],val_y[0:288,:]\n",
        "            print(\"train_X.shape, train_y.shape, test_X.shape, test_y.shape=\",train_X.shape, train_y.shape, test_X.shape, test_y.shape)\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\n",
        "            n_features = train_X.shape[2]\n",
        "            err,err2 = MODEL_EN_DC(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv)\n",
        "            allres = np.array([item[0] for item in err])\n",
        "            allres2 = np.array([item[0] for item in err2])\n",
        "            for item in allres:\n",
        "              item = round(item,4)\n",
        "              outfile_en_dc.write(str(item) + \"  \")\n",
        "            for item in allres2:\n",
        "              item = round(item,4)\n",
        "              train_en_dc.write(str(item) + \"  \")\n",
        "            #np.savetxt(outfile_en_dc,  allres  , fmt='%1.4f', newline=' ' )\n",
        "            #np.savetxt(train_en_dc,  allres2  , fmt='%1.4f', newline=' ' )\n",
        "        outfile_en_dc.write(\"\\n\")    \n",
        "        train_en_dc.write(\"\\n\") \n",
        "            \n",
        "outfile_en_dc.close()    \n",
        "train_en_dc.close()        \n",
        "time2 = time.time()\n",
        "print(\"MODEL_ENC_DCD:\",time2-time1)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kEBN_kGRe1Y-"
      },
      "source": [
        "Plotting different bars:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-kpHmUEYb7F8"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import datetime\n",
        "def show_results(r,method):\n",
        "    x = datetime.datetime.now()\n",
        "    print(\"r = \", r)\n",
        "    u=np.mean(r,axis=0)\n",
        "    sigma=np.std(r,axis=0)\n",
        "    lower_bar= u-1.96*(sigma/(np.sqrt(1)))\n",
        "    upper_bar= u+1.96*(sigma/(np.sqrt(1)))\n",
        "    barWidth = 0.3\n",
        "    bars1=u\n",
        "    rl= np.arange(len(bars1))+1\n",
        "    plt.bar(rl, bars1, width = barWidth, color = 'red', edgecolor = 'black', yerr=1.96*(sigma/(np.sqrt(r.shape[0]))), capsize=7, label=method+' with using 10 multivariate steps for 30 iterations')\n",
        "    plt.ylabel('mean of normalized RMSE for 30 iterations')\n",
        "    plt.xlabel('Number of steps')\n",
        "    plt.legend()\n",
        "    plt.show()\n",
        "    "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z0hjyfWgb7F9"
      },
      "source": [
        "\n",
        "ladam = []\n",
        "# change file directory for different stations accordingly\n",
        "file1 = open('results/AnandVihar/results_adam.txt','r')\n",
        "Lines = file1.readlines()\n",
        "for line in Lines:\n",
        "  strline = line.strip()\n",
        "  li = strline.split()\n",
        "  li = li[:10]\n",
        "  ladam.append([float(li[k]) for k in range(10)])\n",
        "\n",
        "dADAM=np.array(ladam)\n",
        "print(dADAM.shape)\n",
        "show_results(dADAM,'fnn_ADAM')\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Mwcp038lb7F-"
      },
      "source": [
        "\n",
        "dLSTM = []\n",
        "file1 = open('results/AnandVihar/results_lstm.txt','r')\n",
        "Lines = file1.readlines()\n",
        "for line in Lines:\n",
        "  strline = line.strip()\n",
        "  li = strline.split()\n",
        "  li = li[:10]\n",
        "  dLSTM.append([float(li[k]) for k in range(10)])\n",
        "  \n",
        "dLSTM=np.array(dLSTM)\n",
        "print(dLSTM.shape)\n",
        "\n",
        "show_results(dLSTM,'LSTM')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EyqiKuBxb7F_"
      },
      "source": [
        "dEN = []\n",
        "file1 = open('results/AnandVihar/results_en_dc.txt','r')\n",
        "Lines = file1.readlines()\n",
        "for line in Lines:\n",
        "  strline = line.strip()\n",
        "  li = strline.split()\n",
        "  li = li[:10]\n",
        "  dEN.append([float(li[k]) for k in range(10)])\n",
        "dEN=np.array(dEN)\n",
        "print(dEN.shape)\n",
        "show_results(dEN,'Encoder-Decoder')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0lmnDaaub7GA"
      },
      "source": [
        "dBI = []\n",
        "file1 = open('results/AnandVihar/results_bi_lstm.txt','r')\n",
        "Lines = file1.readlines()\n",
        "for line in Lines:\n",
        "  strline = line.strip()\n",
        "  li = strline.split()\n",
        "  li = li[:10]\n",
        "  dBI.append([float(li[k]) for k in range(10)])\n",
        "dBI=np.array(dBI)\n",
        "print(dBI.shape)\n",
        "show_results(dBI,'Bi directional LSTM')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rVyqOSA7b7GA"
      },
      "source": [
        "def show_results_Combined(r1,r2,r3,r4):\n",
        "    x = datetime.datetime.now()\n",
        "    print(r1.shape)\n",
        "    u1=np.mean(r1,axis=0)\n",
        "    u2=np.mean(r2,axis=0)\n",
        "    u3=np.mean(r3,axis=0)\n",
        "    u4=np.mean(r4,axis=0)\n",
        "    sigma1=np.std(r1,axis=0)\n",
        "    sigma2=np.std(r2,axis=0)\n",
        "    sigma3=np.std(r3,axis=0)\n",
        "    sigma4=np.std(r4,axis=0)\n",
        "    barWidth = 0.17\n",
        "    bars1=u1\n",
        "    bars2=u2\n",
        "    bars3=u3\n",
        "    bars4=u4\n",
        "    ra1= np.arange(len(bars1))+1\n",
        "    ra2 = [y + barWidth for y in ra1]\n",
        "    ra3 = [y + barWidth for y in ra2]\n",
        "    ra4 = [y + barWidth for y in ra3]\n",
        "    plt.bar(ra1, bars1, width = barWidth, color = 'blue', edgecolor = 'black', yerr=1.96*(sigma1/(np.sqrt(r1.shape[0]))), capsize=7, label='FNN-ADAM')\n",
        "    plt.bar(ra2, bars2, width = barWidth, color = 'yellow', edgecolor = 'black', yerr=1.96*(sigma2/(np.sqrt(r2.shape[0]))), capsize=7, label='LSTM')\n",
        "    plt.bar(ra3, bars3, width = barWidth, color = 'green', edgecolor = 'black', yerr=1.96*(sigma3/(np.sqrt(r3.shape[0]))), capsize=7, label='BD-LSTM')\n",
        "    plt.bar(ra4, bars4, width = barWidth, color = 'red', edgecolor = 'black', yerr=1.96*(sigma4/(np.sqrt(r4.shape[0]))), capsize=7, label='ED-LSTM')\n",
        "    plt.ylabel('mean of RMSE for 30 iterations')\n",
        "    plt.xlabel('Number of steps')\n",
        "    plt.legend()\n",
        "    plt.show()\n",
        "show_results_Combined(dADAM,dLSTM,dBI,dEN)\n",
        "## Change the fontsize of legends according to convenience"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "I_2szcMnb7GB",
        "scrolled": true
      },
      "source": [
        "import plotly.graph_objects as go\n",
        "\n",
        "u1=np.mean(dADAM,axis=0)\n",
        "u2=np.mean(dLSTM,axis=0)\n",
        "u3=np.mean(dEN,axis=0)\n",
        "u4=np.mean(dBI,axis=0)\n",
        "sigma1=np.array([1.96*(np.std(dADAM[:,k])/(np.sqrt(dADAM[:,k].shape[0]))) for k in range(10)])\n",
        "sigma2=np.array([1.96*(np.std(dLSTM[:,k])/(np.sqrt(dLSTM[:,k].shape[0]))) for k in range(10)])\n",
        "sigma3=np.array([1.96*(np.std(dEN[:,k])/(np.sqrt(dEN[:,k].shape[0]))) for k in range(10)])\n",
        "sigma4=np.array([1.96*(np.std(dBI[:,k])/(np.sqrt(dBI[:,k].shape[0]))) for k in range(10)])\n",
        "fig = go.Figure()\n",
        "df=pd.DataFrame()\n",
        "df['ADAM']=u1\n",
        "df['LSTM']=u2\n",
        "df['EN']=u3\n",
        "df['BI']=u4\n",
        "df['steps']=[1,2,3,4,5,6,7,8,9,10]\n",
        "fig.add_trace(go.Bar(x=df['steps'],\n",
        "                y=df['ADAM'],\n",
        "                     error_y=dict(type='data', array=sigma1),\n",
        "                name='FNN-ADAM',\n",
        "                marker_color='rgb(55, 83, 109)'\n",
        "                ))\n",
        "fig.add_trace(go.Bar(x=df['steps'],\n",
        "                y=df['LSTM'],\n",
        "                error_y=dict(type='data', array=sigma2),\n",
        "                name='LSTM',\n",
        "                marker_color='rgb(200,125,150)'\n",
        "                ))\n",
        "fig.add_trace(go.Bar(x=df['steps'],\n",
        "                y=df['EN'],\n",
        "                error_y=dict(type='data', array=sigma3),\n",
        "                name='ED-LSTM',\n",
        "                marker_color='rgb(13, 255, 20)'\n",
        "                ))\n",
        "fig.add_trace(go.Bar(x=df['steps'],\n",
        "                y=df['BI'],\n",
        "                error_y=dict(type='data', array=sigma4),\n",
        "                name='BD-LSTM',\n",
        "                marker_color='rgb(255, 10, 10)'\n",
        "                ))\n",
        "\n",
        "fig.update_layout(\n",
        "    #title='mean of RMSE for 30 iterations' ,\n",
        "    width = 800,\n",
        "    xaxis_tickfont_size=20,\n",
        "    xaxis=dict(\n",
        "        title='Step Size',\n",
        "        titlefont_size=20,\n",
        "        tickfont_size=20,\n",
        "    ),\n",
        "    yaxis=dict(\n",
        "        title='RMSE(Mean)',\n",
        "        titlefont_size=20,\n",
        "        tickfont_size=20,\n",
        "    ),\n",
        "    legend=dict(\n",
        "        x=0,\n",
        "        y=1.0,\n",
        "        font=dict(\n",
        "            size=18\n",
        "        ),\n",
        "        bgcolor='rgba(255, 255, 255, 0)',\n",
        "        bordercolor='rgba(255, 255, 255, 0)'\n",
        "    ),\n",
        "    barmode='group',\n",
        "    bargap=0.25, # gap between bars of adjacent location coordinates.\n",
        "    bargroupgap=0.05# gap between bars of the same location coordinate.\n",
        "   \n",
        "    \n",
        ")\n",
        "fig.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hcxoWaKygYLd"
      },
      "source": [
        "Data Preparation 1: Random Shuffling"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QGwpe6RzgQ-N"
      },
      "source": [
        "\r\n",
        "def Shuffled_BiLSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv):\r\n",
        "    model = Sequential()\r\n",
        "    model.add(Bidirectional(LSTM(50, activation='relu', input_shape=(n_hours,n_features))))\r\n",
        "    model.add(Dense(n_out))\r\n",
        "    model.compile(optimizer='adam', loss='mse')\r\n",
        "# fit model\r\n",
        "    model.fit(train_X, train_y, epochs=64, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle= True)\r\n",
        "    yhat11 = model.predict(test_X, verbose=0)\r\n",
        "    yhat11_inv = ((maxv - minv) * yhat11  + minv)\r\n",
        "    test_y_inv = ((maxv - minv) * test_y  + minv)\r\n",
        "    print(\"Predicted = \", yhat11_inv)\r\n",
        "    print(\"Actual = \", test_y_inv)\r\n",
        "    err = np.sqrt(((np.sum((yhat11_inv-test_y_inv)**2,axis=0))))\r\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step is =\",err )\r\n",
        "    yhat12 = model.predict(train_X, verbose = 0)\r\n",
        "    yhat12_inv = ((maxv - minv) * yhat12  + minv)\r\n",
        "    train_y_inv = ((maxv - minv) * train_y  + minv)\r\n",
        "    err2 = np.sqrt(((np.sum((yhat12_inv-train_y_inv)**2,axis=0))))\r\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step for train is =\",err2 )\r\n",
        "    return err , err2\r\n",
        "\r\n",
        "\r\n",
        "time1 = time.time()\r\n",
        "outfile_bilstm_shuffled = open('results/AnandVihar/results_bilstm_shuffled.txt','a+')\r\n",
        "train_shuffled_bilstm = open('results/AnandVihar/bilstm_shuffled_train.txt','a+')\r\n",
        "for i,k in enumerate(out_range):\r\n",
        "    for j in range(1,31):\r\n",
        "        for m in range(5,6):\r\n",
        "            n_out=k\r\n",
        "            n_hours=m\r\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\r\n",
        "            \r\n",
        "            train_X, train_y = split_sequences(scaled_train,n_hours,n_out)\r\n",
        "            val_X, val_y = split_sequences(scaled_test,n_hours,n_out)\r\n",
        "            print(val_X.shape,val_y.shape)\r\n",
        "            test_X, test_y = val_X[288:,:,:],val_y[288:,:]\r\n",
        "            val_X, val_y = val_X[0:288,:,:],val_y[0:288,:]\r\n",
        "            print(\"train_X.shape, train_y.shape, test_X.shape, test_y.shape=\",train_X.shape, train_y.shape, test_X.shape, test_y.shape)\r\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\r\n",
        "            n_features = train_X.shape[2]\r\n",
        "            print(\"n_features = \",n_features)\r\n",
        "            err , err2= Shuffled_BiLSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv)\r\n",
        "            allres = np.array(err)\r\n",
        "            allres2 = np.array(err2)\r\n",
        "            for item in allres:\r\n",
        "              item = round(item,4)\r\n",
        "              outfile_bilstm_shuffled.write(str(item) + \"  \")\r\n",
        "            for item in allres2:\r\n",
        "              item = round(item,4)\r\n",
        "              train_shuffled_bilstm.write(str(item) + \"  \")\r\n",
        "            #np.savetxt(outfile_bilstm_shuffled,  allres  , fmt='%1.4f', newline=' ' )\r\n",
        "            #np.savetxt(train_shuffled_bilstm,  allres2  , fmt='%1.4f', newline=' ' )\r\n",
        "        outfile_bilstm_shuffled.write(\"\\n\")    \r\n",
        "        train_shuffled_bilstm.write(\"\\n\")     \r\n",
        "\r\n",
        "outfile_bilstm_shuffled.close()    \r\n",
        "train_shuffled_bilstm.close()            \r\n",
        "time2 = time.time()        \r\n",
        "print(\"MODEL_time:\",time2-time1)\r\n",
        "#np.savetxt(\"10i_\"+str(result.shape[1])+\"o_FNN_adam_\"+str(result.shape[0])+\"i.csv\"+str(datetime.datetime.now()), result, delimiter=\",\")\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NORv6MwCzYZ4"
      },
      "source": [
        "Data Preparation 2: Univariate vs Multivariate Analysis using our best model- Bi-LSTM"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RZa8No3qzXy4"
      },
      "source": [
        "\r\n",
        "def split_sequences_univariate(sequences, n_steps_in, n_steps_out):\r\n",
        "    X, y = list(), list()\r\n",
        "    for i in range(len(sequences)):\r\n",
        "    # find the end of this pattern\r\n",
        "        end_ix = i + n_steps_in\r\n",
        "        out_end_ix = end_ix + n_steps_out-1\r\n",
        "    # check if we are beyond the dataset\r\n",
        "        if out_end_ix > len(sequences):\r\n",
        "            break\r\n",
        "        # gather input and output parts of the pattern\r\n",
        "        \r\n",
        "        seq_x, seq_y = sequences[i:end_ix], sequences[end_ix-1:out_end_ix]\r\n",
        "        X.append(seq_x)\r\n",
        "        y.append(seq_y)\r\n",
        "    return array(X), array(y)\r\n",
        "\r\n",
        "def MODEL_Univariate_Bi_LSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv):\r\n",
        "    model = Sequential()\r\n",
        "    model.add(Bidirectional(LSTM(50, activation='relu', input_shape=(n_hours,n_features))))\r\n",
        "    model.add(Dense(n_out))\r\n",
        "    model.compile(optimizer='adam', loss='mse')\r\n",
        "    model.fit(train_X, train_y, epochs=1000, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle=False)\r\n",
        "    yhat11 = model.predict(test_X, verbose=0)\r\n",
        "    yhat11_inv = ((maxv - minv) * yhat11  + minv)\r\n",
        "    test_y_inv = ((maxv - minv) * test_y  + minv)\r\n",
        "    print(\"Predicted = \", yhat11_inv)\r\n",
        "    print(\"Actual = \", test_y_inv)\r\n",
        "    err = np.sqrt(((np.sum((yhat11_inv-test_y_inv)**2,axis=0))))\r\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step is =\",err )\r\n",
        "    yhat12 = model.predict(train_X, verbose = 0)\r\n",
        "    yhat12_inv = ((maxv - minv) * yhat12  + minv)\r\n",
        "    train_y_inv = ((maxv - minv) * train_y  + minv)\r\n",
        "    err2 = np.sqrt(((np.sum((yhat12_inv-train_y_inv)**2,axis=0))))\r\n",
        "    print(\"RMSE error using \",n_hours,\" state size to predict \",n_out,\" step for train is =\",err2 )\r\n",
        "    return err , err2\r\n",
        "\r\n",
        "\r\n",
        "\r\n",
        "dataset_train = read_csv('train_pollution.csv', header=0, index_col=0)\r\n",
        "dataset_test = read_csv('test_pollution.csv', header=0, index_col=0)\r\n",
        "values_train = dataset_train.values\r\n",
        "values_test = dataset_test.values\r\n",
        "values_train=values_train[:,-7]\r\n",
        "values_test=values_test[:,-7]\r\n",
        "values_train = values_train.astype('float32')\r\n",
        "values_test = values_test.astype('float32')\r\n",
        "scaler = MinMaxScaler(feature_range=(0, 1))\r\n",
        "concat = np.concatenate((values_train,values_test),axis = 0)\r\n",
        "print(concat.shape)\r\n",
        "maxv = max(concat)\r\n",
        "minv = min(concat)\r\n",
        "scaled = (concat-minv)/(maxv-minv)\r\n",
        "scaled_train = scaled[:values_train.shape[0]]\r\n",
        "scaled_test = scaled[values_train.shape[0]:]\r\n",
        "outfile_bi_lstm = open('results/AnandVihar/results_bi_lstm_univariate.txt','a+')\r\n",
        "train_bi_lstm = open('results/AnandVihar/bi_lstm_univarate_train.txt','a+')\r\n",
        "time1 = time.time()\r\n",
        "for i,k in enumerate(out_range):\r\n",
        "    for j in range(1,31):\r\n",
        "        for m in range(5,6):\r\n",
        "            n_out=k\r\n",
        "            n_hours=m\r\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\r\n",
        "            \r\n",
        "            train_X, train_y = split_sequences_univariate(scaled_train,n_hours,n_out)\r\n",
        "            val_X, val_y = split_sequences_univariate(scaled_test,n_hours,n_out)\r\n",
        "            val_X = val_X.reshape(val_X.shape[0] , val_X.shape[1] , 1)\r\n",
        "            train_X = train_X.reshape(train_X.shape[0] , train_X.shape[1] , 1)\r\n",
        "            print(val_X.shape,val_y.shape)\r\n",
        "            test_X, test_y = val_X[288:,:,:],val_y[288:,:]\r\n",
        "            val_X, val_y = val_X[0:288,:,:],val_y[0:288,:]\r\n",
        "            print(\"train_X.shape, train_y.shape, test_X.shape, test_y.shape=\",train_X.shape, train_y.shape, test_X.shape, test_y.shape)\r\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\r\n",
        "            n_features = 1\r\n",
        "            err , err2=MODEL_Univariate_Bi_LSTM(train_X,test_X,train_y,test_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv)\r\n",
        "            allres = np.array(err)\r\n",
        "            allres2 = np.array(err2)\r\n",
        "            for item in allres:\r\n",
        "              item = round(item,4)\r\n",
        "              outfile_bi_lstm.write(str(item) + \"  \")\r\n",
        "            for item in allres2:\r\n",
        "              item = round(item,4)\r\n",
        "              train_bi_lstm.write(str(item) + \"  \")\r\n",
        "            #np.savetxt(outfile_bi_lstm,  allres  , fmt='%1.4f', newline=' ' )\r\n",
        "            #np.savetxt(train_bi_lstm,  allres2  , fmt='%1.4f', newline=' ' )\r\n",
        "        outfile_bi_lstm.write(\"\\n\")    \r\n",
        "        train_bi_lstm.write(\"\\n\") \r\n",
        "\r\n",
        "outfile_bi_lstm.close()    \r\n",
        "train_bi_lstm.close()            \r\n",
        "        \r\n",
        "time2 = time.time()\r\n",
        "print(\"MODEL_Bi_LSTM:\",time2-time1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mruxaASuADaO"
      },
      "source": [
        "shuffled = []\r\n",
        "file1 = open('results/AnandVihar/results_bilstm_shuffled.txt','r')\r\n",
        "#file1 = open('results/AnandVihar/1_shuffled_test.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  shuffled.append([float(li[k]) for k in range(10)])\r\n",
        "shuffled = np.array(shuffled)\r\n",
        "flat_shuffled = np.ravel(shuffled , order = 'F')\r\n",
        "print(\"Mean = \", np.mean(flat_shuffled))\r\n",
        "print(\"Error = \", 1.96*(np.std(flat_shuffled)/(np.sqrt(flat_shuffled.shape[0]))))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p-8jRM5SCF6Z"
      },
      "source": [
        "Plotting for results of different models"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lxAsp7OXbuaw"
      },
      "source": [
        "##Get the test results\r\n",
        "#FNN-ADAM\r\n",
        "dADAM_train = []\r\n",
        "file1 = open('results/AnandVihar/adam_train.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  dADAM_train.append([float(li[k]) for k in range(10)])\r\n",
        "dADAM_train=np.array(dADAM_train)\r\n",
        "\r\n",
        "#LSTM\r\n",
        "dLSTM_train = []\r\n",
        "file1 = open('results/AnandVihar/lstm_train.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  dLSTM_train.append([float(li[k]) for k in range(10)])\r\n",
        "dLSTM_train=np.array(dLSTM_train)\r\n",
        "\r\n",
        "#ED-LSTM\r\n",
        "dEN_train = []\r\n",
        "file1 = open('results/AnandVihar/en_dc_train.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  dEN_train.append([float(li[k]) for k in range(10)])\r\n",
        "dEN_train=np.array(dEN_train)\r\n",
        "\r\n",
        "#BD-LSTM\r\n",
        "dBI_train = []\r\n",
        "file1 = open('results/AnandVihar/bi_lstm_train.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  dBI_train.append([float(li[k]) for k in range(10)])\r\n",
        "dBI_train=np.array(dBI_train)\r\n",
        "\r\n",
        "#BD-LSTM*\r\n",
        "shuffled_train = []\r\n",
        "file1 = open('results/AnandVihar/bilstm_shuffled_train.txt','r')\r\n",
        "#file1 = open('results/AnandVihar/1_shuffled_train.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  shuffled_train.append([float(li[k]) for k in range(10)])\r\n",
        "shuffled_train=np.array(shuffled_train)\r\n",
        "\r\n",
        "#UBD-LSTM\r\n",
        "UBD_train = []\r\n",
        "file1 = open('results/AnandVihar/bi_lstm_univarate_train.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  UBD_train.append([float(li[k]) for k in range(10)])\r\n",
        "UBD_train=np.array(UBD_train)\r\n",
        "#SBD-LSTM\r\n",
        "SBD_train = []\r\n",
        "file1 = open('results/AnandVihar/seasonal_bi_lstm_train.txt','r')\r\n",
        "#file1 = open('results/AnandVihar/1_seasonal_train.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  SBD_train.append([float(li[k]) for k in range(10)])\r\n",
        "SBD_train=np.array(SBD_train)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VYyUxTXDCEcz"
      },
      "source": [
        "## Corresponding Test RMSE with Mean and Error across 10 horizons\r\n",
        "#flat_dSGD = np.ravel(dSGD , order = 'F')\r\n",
        "flat_dADAM = np.ravel(dADAM , order = 'F')\r\n",
        "flat_dLSTM = np.ravel(dLSTM , order = 'F')\r\n",
        "flat_dEN = np.ravel(dEN , order = 'F')\r\n",
        "flat_dBI = np.ravel(dBI , order = 'F')\r\n",
        "\r\n",
        "dBI_univ = []\r\n",
        "file1 = open('results/AnandVihar/results_bi_lstm_univariate.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  dBI_univ.append([float(li[k]) for k in range(10)])\r\n",
        "dBI_univ = np.array(dBI_univ)\r\n",
        "flat_dBI_univ = np.ravel(dBI_univ , order = 'F')\r\n",
        "\r\n",
        "dBI_season = []\r\n",
        "file1 = open('results/AnandVihar/result_seasonal_bi_lstm.txt','r')\r\n",
        "#file1 = open('results/AnandVihar/1_seasonal_test.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:10]\r\n",
        "  dBI_season.append([float(li[k]) for k in range(10)])\r\n",
        "dBI_season = np.array(dBI_season)\r\n",
        "flat_dBI_season = np.ravel(dBI_season , order = 'F')\r\n",
        "\r\n",
        "#m1 = np.mean(flat_dSGD)\r\n",
        "m2 = np.mean(flat_dADAM)\r\n",
        "m3 = np.mean(flat_dLSTM)\r\n",
        "m4 = np.mean(flat_dEN)\r\n",
        "m5 = np.mean(flat_dBI)\r\n",
        "m6 = np.mean(flat_dBI_univ)\r\n",
        "m7 = np.mean(flat_shuffled)\r\n",
        "m8 = np.mean(flat_dBI_season)\r\n",
        "\r\n",
        "#s1 = 1.96*(np.std(flat_dSGD)/(np.sqrt(flat_dSGD.shape[0])))\r\n",
        "s2 = 1.96*(np.std(flat_dADAM)/(np.sqrt(flat_dADAM.shape[0])))\r\n",
        "s3 = 1.96*(np.std(flat_dLSTM)/(np.sqrt(flat_dLSTM.shape[0])))\r\n",
        "s4 = 1.96*(np.std(flat_dEN)/(np.sqrt(flat_dEN.shape[0])))\r\n",
        "s5 = 1.96*(np.std(flat_dBI)/(np.sqrt(flat_dBI.shape[0])))\r\n",
        "s6 = 1.96*(np.std(flat_dBI_univ)/(np.sqrt(flat_dBI_univ.shape[0])))\r\n",
        "s7 = 1.96*(np.std(flat_shuffled)/(np.sqrt(flat_shuffled.shape[0])))\r\n",
        "s8 = 1.96*(np.std(flat_dBI_season)/(np.sqrt(flat_dBI_season.shape[0])))\r\n",
        "\r\n",
        "## Corresponding Train RMSE with Mean and Error across 10 horizons\r\n",
        "flat_dADAM_train = np.ravel(dADAM_train , order = 'F')\r\n",
        "flat_dLSTM_train = np.ravel(dLSTM_train , order = 'F')\r\n",
        "flat_dEN_train = np.ravel(dEN_train , order = 'F')\r\n",
        "flat_dBI_train = np.ravel(dBI_train , order = 'F')\r\n",
        "flat_shuffled_train = np.ravel(shuffled_train , order = 'F')\r\n",
        "flat_UBD_train = np.ravel(UBD_train , order = 'F')\r\n",
        "flat_SBD_train = np.ravel(SBD_train , order = 'F')\r\n",
        "\r\n",
        "n2 = np.mean(flat_dADAM_train)\r\n",
        "n3 = np.mean(flat_dLSTM_train)\r\n",
        "n4 = np.mean(flat_dEN_train)\r\n",
        "n5 = np.mean(flat_dBI_train)\r\n",
        "n6 = np.mean(flat_UBD_train)\r\n",
        "n7 = np.mean(flat_shuffled_train)\r\n",
        "n8 = np.mean(flat_SBD_train)\r\n",
        "\r\n",
        "t2 = 1.96*(np.std(flat_dADAM_train)/(np.sqrt(flat_dADAM_train.shape[0])))\r\n",
        "t3 = 1.96*(np.std(flat_dLSTM_train)/(np.sqrt(flat_dLSTM_train.shape[0])))\r\n",
        "t4 = 1.96*(np.std(flat_dEN_train)/(np.sqrt(flat_dEN_train.shape[0])))\r\n",
        "t5 = 1.96*(np.std(flat_dBI_train)/(np.sqrt(flat_dBI_train.shape[0])))\r\n",
        "t6 = 1.96*(np.std(flat_UBD_train)/(np.sqrt(flat_UBD_train.shape[0])))\r\n",
        "t7 = 1.96*(np.std(flat_shuffled_train)/(np.sqrt(flat_shuffled_train.shape[0])))\r\n",
        "t8 = 1.96*(np.std(flat_SBD_train)/(np.sqrt(flat_SBD_train.shape[0])))\r\n",
        "\r\n",
        "X = ['Train','Test']\r\n",
        "fig = go.Figure()\r\n",
        "\r\n",
        "fig.add_trace(go.Bar(x= X,\r\n",
        "                y=np.array([n2,m2]),\r\n",
        "                     error_y=dict(type='data', array= np.array([t2,s2])),\r\n",
        "                name='FNN-ADAM',\r\n",
        "                marker_color='rgb(55, 83, 109)'\r\n",
        "                ))\r\n",
        "fig.add_trace(go.Bar(x= X,\r\n",
        "                y= np.array([n3,m3]),\r\n",
        "                error_y=dict(type='data', array= np.array([t3,s3])),\r\n",
        "                name='LSTM',\r\n",
        "                marker_color='rgb(200,125,150)'\r\n",
        "                ))\r\n",
        "fig.add_trace(go.Bar(x=X,\r\n",
        "                y= np.array([n4,m4]),\r\n",
        "                error_y=dict(type='data', array= np.array([t4,s4])),\r\n",
        "                name='ED-LSTM',\r\n",
        "                marker_color='rgb(13, 255, 20)'\r\n",
        "                ))\r\n",
        "fig.add_trace(go.Bar(x=X,\r\n",
        "                y= np.array([n5,m5]),\r\n",
        "                error_y=dict(type='data', array= np.array([t5,s5])),\r\n",
        "                name='BD-LSTM',\r\n",
        "                marker_color='rgb(255, 10, 10)'\r\n",
        "                ))\r\n",
        "\r\n",
        "\r\n",
        "fig.update_layout(\r\n",
        "    #title='Mean and Error of RMSE for 10 prediction horizons' ,\r\n",
        "    xaxis_tickfont_size=20,\r\n",
        "    width = 750,\r\n",
        "    yaxis=dict(\r\n",
        "        title='RMSE(Mean)',\r\n",
        "        titlefont_size=20,\r\n",
        "        tickfont_size=20,\r\n",
        "    ),\r\n",
        "    legend=dict(\r\n",
        "        x=0.94,\r\n",
        "        y=1.0,\r\n",
        "        font = dict(size = 18),\r\n",
        "        bgcolor='rgba(255, 255, 255, 0)',\r\n",
        "        bordercolor='rgba(255, 255, 255, 0)'\r\n",
        "    ),\r\n",
        "    barmode='group',\r\n",
        "    bargap=0.2, # gap between bars of adjacent location coordinates.\r\n",
        "    bargroupgap=0.05# gap between bars of the same location coordinate.\r\n",
        "   \r\n",
        ")\r\n",
        "fig.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qEWVoKoKFyyc"
      },
      "source": [
        "X = ['Train','Test']\r\n",
        "fig = go.Figure()\r\n",
        "fig.add_trace(go.Bar(x=X,\r\n",
        "                y= np.array([n5,m5]),\r\n",
        "                error_y=dict(type='data', array= np.array([t5,s5])),\r\n",
        "                name='BD-LSTM',\r\n",
        "                marker_color='rgb(255, 10, 10)'\r\n",
        "                ))\r\n",
        "fig.add_trace(go.Bar(x= X,\r\n",
        "                y=np.array([n7,m7]),\r\n",
        "                     error_y=dict(type='data', array= np.array([t7,s7])),\r\n",
        "                name='BD-LSTM*',\r\n",
        "                marker_color='#FF851B'\r\n",
        "                ))\r\n",
        "fig.add_trace(go.Bar(x=X,\r\n",
        "                y= np.array([n6,m6]),\r\n",
        "                error_y=dict(type='data', array= np.array([t6,s6])),\r\n",
        "                name='UBD-LSTM',\r\n",
        "                marker_color='lightseagreen'\r\n",
        "                ))\r\n",
        "            \r\n",
        "fig.add_trace(go.Bar(x=X,\r\n",
        "                y= np.array([n8,m8]),\r\n",
        "                error_y=dict(type='data', array= np.array([t8,s8])),\r\n",
        "                name='SBD-LSTM',\r\n",
        "                marker_color='darkblue'\r\n",
        "                ))\r\n",
        "\r\n",
        "\r\n",
        "fig.update_layout(\r\n",
        "    #title='Mean and Error of RMSE for 10 prediction horizons' ,\r\n",
        "    xaxis_tickfont_size=20,\r\n",
        "    width = 700,\r\n",
        "    yaxis=dict(\r\n",
        "        title='RMSE(Mean)',\r\n",
        "        titlefont_size=20,\r\n",
        "        tickfont_size=20,\r\n",
        "    ),\r\n",
        "    legend=dict(\r\n",
        "        x=0.94,\r\n",
        "        y=1.0,\r\n",
        "        font = dict(size = 18),\r\n",
        "        bgcolor='rgba(255, 255, 255, 0)',\r\n",
        "        bordercolor='rgba(255, 255, 255, 0)'\r\n",
        "    ),\r\n",
        "    barmode='group',\r\n",
        "    bargap=0.2, # gap between bars of adjacent location coordinates.\r\n",
        "    bargroupgap=0.05# gap between bars of the same location coordinate.\r\n",
        "   \r\n",
        ")\r\n",
        "fig.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "92PqLLeGi-AV"
      },
      "source": [
        "Table for all time steps\r\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "s8vOr2AfjIo4"
      },
      "source": [
        "row_train_mean = [n2,n3,n4,n5,n7,n6,n8,t2,t3,t4,t5,t7,t6,t8]\r\n",
        "row_test_mean = [m2,m3,m4,m5,m7,m6,m8,s2,s3,s4,s5,s7,s6,s8]\r\n",
        "all_meth = [dADAM, dLSTM, dEN, dBI, shuffled, dBI_univ, dBI_season]\r\n",
        "step1_mean = [np.mean(arr[:,0]) for arr in all_meth]\r\n",
        "step1_err = [1.96*(np.std(arr[:,0])/(np.sqrt(arr[:,0].shape[0]))) for arr in all_meth]\r\n",
        "step2_mean = [np.mean(arr[:,1]) for arr in all_meth]\r\n",
        "step2_err = [1.96*(np.std(arr[:,1])/(np.sqrt(arr[:,1].shape[0]))) for arr in all_meth]\r\n",
        "step3_mean = [np.mean(arr[:,2]) for arr in all_meth]\r\n",
        "step3_err = [1.96*(np.std(arr[:,2])/(np.sqrt(arr[:,2].shape[0]))) for arr in all_meth]\r\n",
        "step4_mean = [np.mean(arr[:,3]) for arr in all_meth]\r\n",
        "step4_err = [1.96*(np.std(arr[:,3])/(np.sqrt(arr[:,3].shape[0]))) for arr in all_meth]\r\n",
        "step5_mean = [np.mean(arr[:,4]) for arr in all_meth]\r\n",
        "step5_err = [1.96*(np.std(arr[:,4])/(np.sqrt(arr[:,4].shape[0]))) for arr in all_meth]\r\n",
        "step6_mean = [np.mean(arr[:,5]) for arr in all_meth]\r\n",
        "step6_err = [1.96*(np.std(arr[:,5])/(np.sqrt(arr[:,5].shape[0]))) for arr in all_meth]\r\n",
        "step7_mean = [np.mean(arr[:,6]) for arr in all_meth]\r\n",
        "step7_err = [1.96*(np.std(arr[:,6])/(np.sqrt(arr[:,6].shape[0]))) for arr in all_meth]\r\n",
        "step8_mean = [np.mean(arr[:,7]) for arr in all_meth]\r\n",
        "step8_err = [1.96*(np.std(arr[:,7])/(np.sqrt(arr[:,7].shape[0]))) for arr in all_meth]\r\n",
        "step9_mean = [np.mean(arr[:,8]) for arr in all_meth]\r\n",
        "step9_err = [1.96*(np.std(arr[:,8])/(np.sqrt(arr[:,8].shape[0]))) for arr in all_meth]\r\n",
        "step10_mean = [np.mean(arr[:,9]) for arr in all_meth]\r\n",
        "step10_err = [1.96*(np.std(arr[:,9])/(np.sqrt(arr[:,9].shape[0]))) for arr in all_meth]\r\n",
        "##print the corresponding list to get mean/error"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ijjsdtMcA4xQ"
      },
      "source": [
        "1 months ahead prediction using Univariate UBD-LSTM model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qcQwl7kvgHET"
      },
      "source": [
        "\r\n",
        "def split_sequences_univariate(sequences, n_steps_in, n_steps_out):\r\n",
        "    X, y = list(), list()\r\n",
        "    for i in range(len(sequences)):\r\n",
        "    # find the end of this pattern\r\n",
        "        end_ix = i + n_steps_in\r\n",
        "        out_end_ix = end_ix + n_steps_out-1\r\n",
        "    # check if we are beyond the dataset\r\n",
        "        if out_end_ix > len(sequences):\r\n",
        "            break\r\n",
        "        # gather input and output parts of the pattern\r\n",
        "        \r\n",
        "        seq_x, seq_y = sequences[i:end_ix], sequences[end_ix-1:out_end_ix]\r\n",
        "        X.append(seq_x)\r\n",
        "        y.append(seq_y)\r\n",
        "    return array(X), array(y)\r\n",
        "\r\n",
        "def monthly(train_X,train_y,val_X,val_y,n_out,n_hours,n_features):\r\n",
        "    model = Sequential()\r\n",
        "    model.add(Bidirectional(LSTM(50, activation='relu', input_shape=(n_hours,n_features))))\r\n",
        "    model.add(Dense(n_out))\r\n",
        "    model.compile(optimizer='adam', loss='mse')\r\n",
        "    model.fit(train_X, train_y, epochs=1000, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle= False)\r\n",
        "    return model\r\n",
        "\r\n",
        "def predict_model(model_old,n_hours,n_features,n_out):\r\n",
        "  new_model = Sequential()\r\n",
        "  new_model.add(Bidirectional(LSTM(50, activation='relu', input_shape=(n_hours,n_features))))\r\n",
        "  new_model.add(Dense(n_out))\r\n",
        "  new_model.build(input_shape = (1,n_hours,n_features))\r\n",
        "  old_weights = model_old.get_weights()\r\n",
        "  new_model.set_weights(old_weights)\r\n",
        "  new_model.compile(optimizer='adam', loss='mse')\r\n",
        "  return new_model\r\n",
        "\r\n",
        "dataset_train = read_csv('train_pollution.csv', header=0, index_col=0)\r\n",
        "dataset_test = read_csv('test_pollution.csv', header=0, index_col=0)\r\n",
        "values_train = dataset_train.values\r\n",
        "values_test = dataset_test.values\r\n",
        "values_train=values_train[:,-7]\r\n",
        "values_test=values_test[:,-7]\r\n",
        "values_train = values_train.astype('float32')\r\n",
        "values_test = values_test.astype('float32')\r\n",
        "scaler = MinMaxScaler(feature_range=(0, 1))\r\n",
        "concat = np.concatenate((values_train,values_test),axis = 0)\r\n",
        "print(concat.shape)\r\n",
        "maxv = max(concat)\r\n",
        "minv = min(concat)\r\n",
        "scaled = (concat-minv)/(maxv-minv)\r\n",
        "scaled_train = scaled[:values_train.shape[0]]\r\n",
        "scaled_test = scaled[values_train.shape[0]:]\r\n",
        "time1 = time.time()\r\n",
        "outfile_monthly = open('results/AnandVihar/monthly.txt','a+')\r\n",
        "for i,k in enumerate(out_range):\r\n",
        "    for j in range(1,31):\r\n",
        "        for m in range(5,6):\r\n",
        "            n_out=k\r\n",
        "            n_hours=m\r\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\r\n",
        "            train_X, train_y = split_sequences_univariate(scaled_train,n_hours,n_out)\r\n",
        "            val_X, val_y = split_sequences_univariate(scaled_test,n_hours,n_out)\r\n",
        "            val_X = val_X.reshape(val_X.shape[0] , val_X.shape[1] , 1)\r\n",
        "            train_X = train_X.reshape(train_X.shape[0] , train_X.shape[1] , 1)\r\n",
        "            test_X, test_y = val_X[288:,:,:],val_y[288:,:]\r\n",
        "            val_X, val_y = val_X[0:288,:,:],val_y[0:288,:]\r\n",
        "            print(\"train_X.shape, train_y.shape, test_X.shape, test_y.shape=\",train_X.shape, train_y.shape, test_X.shape, test_y.shape)\r\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\r\n",
        "            n_features = 1\r\n",
        "            print(\"n_features = \",n_features)\r\n",
        "            model_monthly= monthly(train_X,train_y,val_X,val_y,n_out,n_hours,n_features)\r\n",
        "            values_list = []\r\n",
        "            #creating the initial input from last 5 values of December which was distributed to test data\r\n",
        "            #now this last 5 window will be used to create next 30 days of predictions (= 30*3 = 90 values) , 3 values for each day.\r\n",
        "            ip = test_X[-1,:,:]\r\n",
        "            initial_input = ip.reshape(1,ip.shape[0],1)\r\n",
        "            print(\"initial_input=\",initial_input.shape)\r\n",
        "            new_model = predict_model(model_monthly,n_hours,n_features,n_out)\r\n",
        "            for st in range(9):\r\n",
        "              next_10 = new_model.predict(initial_input)\r\n",
        "              values_list.extend(next_10)\r\n",
        "              new_ip = next_10[:,5:]\r\n",
        "              initial_input = new_ip.reshape(1,new_ip.shape[1],1)\r\n",
        "\r\n",
        "            write_list = [inner for outer in values_list for inner in outer]\r\n",
        "            final_write_list = [(maxv - minv)*item +minv for item in write_list]\r\n",
        "            for item in final_write_list:\r\n",
        "              item = round(item,4)\r\n",
        "              outfile_monthly.write(str(item) + \"  \")  \r\n",
        "\r\n",
        "        outfile_monthly.write(\"\\n\")    \r\n",
        "\r\n",
        "outfile_monthly.close()         \r\n",
        "time2 = time.time()        \r\n",
        "print(\"MODEL_time:\",time2-time1)\r\n",
        "\r\n",
        "\r\n",
        "\r\n",
        "\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NWHQGq4pkJOT"
      },
      "source": [
        "Plot for one-month ahead with uncertainty"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "C-wxER2XVDuP"
      },
      "source": [
        "\r\n",
        "from datetime import date, datetime, timedelta\r\n",
        "\r\n",
        "def datetime_range(start, end, delta):\r\n",
        "    current = start\r\n",
        "    if not isinstance(delta, timedelta):\r\n",
        "        delta = timedelta(**delta)\r\n",
        "    while current < end:\r\n",
        "        yield current\r\n",
        "        current += delta\r\n",
        "\r\n",
        "\r\n",
        "start = datetime(2020,12,11)\r\n",
        "end = datetime(2021,1,10)\r\n",
        "time_axis = []\r\n",
        "#this unlocks the following interface:\r\n",
        "for dt in datetime_range(start, end, {'days': 0, 'hours':8}):\r\n",
        "  time_axis.append(dt)\r\n",
        "\r\n",
        "future_prediction = []\r\n",
        "file1 = open('results/AnandVihar/monthly.txt','r')\r\n",
        "Lines = file1.readlines()\r\n",
        "for line in Lines:\r\n",
        "  strline = line.strip()\r\n",
        "  li = strline.split()\r\n",
        "  li = li[:90]\r\n",
        "  future_prediction.append([float(li[k]) for k in range(90)])\r\n",
        "\r\n",
        "future_prediction = np.array(future_prediction)\r\n",
        "print(future_prediction.shape)\r\n",
        "pred_values = future_prediction.mean(axis=0)\r\n",
        "error = np.array([1.96*(np.std(future_prediction[:,k])/(np.sqrt(future_prediction[:,k].shape[0]))) for k in range(90)]) \r\n",
        "per_05 = pred_values + error\r\n",
        "per_95 = pred_values - error\r\n",
        "fig = plt.figure()\r\n",
        "params = {'legend.fontsize': 12}\r\n",
        "plt.rcParams.update(params)\r\n",
        "ax = fig.add_subplot(111) \r\n",
        "mu = ax.plot(time_axis, pred_values, color='blue', label = 'Mean')\r\n",
        "s1 = ax.plot(time_axis, per_05, color='green', alpha=0.4,label = 'Error')\r\n",
        "s2 = ax.plot(time_axis, per_95, color='green', alpha=0.4)\r\n",
        "ax.fill_between(time_axis, per_05, per_95, facecolor='g', alpha=0.4)\r\n",
        "ax.set_ylabel('PM2.5 Concentration') \r\n",
        "ax.set_xlabel('Time') \r\n",
        "plt.xticks(rotation = 45)\r\n",
        "# plt.title('Predicted')\r\n",
        "ax.legend()\r\n",
        "#plt.savefig(\"/content/drive/My Drive/covid project/Model/Results/\"+name+\"/\"+Mname+'/future_pred.png',dpi=300) \r\n",
        "plt.show()\r\n",
        "plt.close()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E-NEVGOpjlIW"
      },
      "source": [
        "Trace plot for month of September using seasonal model SBD-LSTM"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "H8MvNscqj8go"
      },
      "source": [
        "\r\n",
        "def trace_model(train_X,july_X,train_y,july_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv):\r\n",
        "    model = Sequential()\r\n",
        "    model.add(Bidirectional(LSTM(50, activation='relu',input_shape=(n_hours,n_features))))\r\n",
        "    model.add(Dense(n_out))\r\n",
        "    model.compile(optimizer='adam', loss='mse')\r\n",
        "# fit model\r\n",
        "    model.fit(train_X, train_y, epochs=50, batch_size=20, validation_data=(val_X, val_y), verbose=0, shuffle=False)\r\n",
        "\r\n",
        "    july_pred = model.predict(july_X,verbose = 0)\r\n",
        "    july_inv = ((maxv - minv) * july_pred  + minv)\r\n",
        "    return july_inv\r\n",
        "\r\n",
        "\r\n",
        "time1 = time.time()\r\n",
        "dataset_train = read_csv('train_pollution.csv', header=0, index_col=0)\r\n",
        "dataset_test = read_csv('test_pollution.csv', header=0, index_col=0)\r\n",
        "values_train = dataset_train.values\r\n",
        "values_test = dataset_test.values\r\n",
        "values_train=values_train[:,1:-6]\r\n",
        "values_test=values_test[:,1:-6]\r\n",
        "values_train = values_train.astype('float32')\r\n",
        "values_test = values_test.astype('float32')\r\n",
        "scaler = MinMaxScaler(feature_range=(0, 1))\r\n",
        "concat = np.concatenate((values_train,values_test),axis = 0)\r\n",
        "maxv = max(concat[:,-1])\r\n",
        "minv = min(concat[:,-1])\r\n",
        "scaled = scaler.fit_transform(concat)\r\n",
        "scaled_train = scaled[:values_train.shape[0],:]\r\n",
        "scaled_test = scaled[values_train.shape[0]:,:]\r\n",
        "allres = []\r\n",
        "for i,k in enumerate(out_range):\r\n",
        "    for j in range(1,31):\r\n",
        "        for m in range(5,6):\r\n",
        "            n_out=k\r\n",
        "            n_hours=m\r\n",
        "            print(\"Number of steps predicted= \",n_out,\"Number of steps Multivariate threads ussed= \",n_hours, \"iteration number= \",j)\r\n",
        "            \r\n",
        "            train_X, train_y = split_sequences(scaled_train,n_hours,n_out)\r\n",
        "            val_X, val_y = split_sequences(scaled_test,n_hours,n_out)\r\n",
        "            july_X, july_y = split_sequences(scaled_july,n_hours,n_out)\r\n",
        "            print(val_X.shape,val_y.shape)\r\n",
        "            print(\"july shape=\",july_X.shape)\r\n",
        "            #test_X, test_y = val_X[365:,:,:],val_y[365:,:]\r\n",
        "            val_X, val_y = val_X[0:365,:,:],val_y[0:365,:]\r\n",
        "            print(\"train_X.shape, train_y.shape =\",train_X.shape, train_y.shape)\r\n",
        "            print(\"val X, val y\",val_X.shape, val_y.shape)\r\n",
        "            n_features = train_X.shape[2]\r\n",
        "            #july_X = test_X\r\n",
        "            #july_y = test_y\r\n",
        "            july_pred = trace_model(train_X,july_X,train_y,july_y,val_X,val_y,n_out,n_hours,n_features,minv,maxv)\r\n",
        "            july_pred = np.array(july_pred)\r\n",
        "            predicted = []\r\n",
        "            for ctr in range(0,july_pred.shape[0],10):\r\n",
        "              predicted.append(july_pred[ctr])\r\n",
        "            predicted = np.array(predicted)  \r\n",
        "            allres.append(predicted.flatten())\r\n",
        "   \r\n",
        "july_y = ((maxv - minv) * july_y  + minv)\r\n",
        "actual = []\r\n",
        "for ctr in range(0,july_y.shape[0],10):\r\n",
        "  actual.append(july_y[ctr])\r\n",
        "actual = np.array(actual).flatten()\r\n",
        "allres = np.array(allres)\r\n",
        "final_pred_y = np.mean(allres, axis = 0)\r\n",
        "\r\n",
        "time2 = time.time()\r\n",
        "print(\"MODEL_july:\",time2-time1)\r\n",
        "\r\n",
        "## Plotting\r\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qpFdMokVkBWv"
      },
      "source": [
        "Corresponding Trace Plot"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LBN4o3ZowxGP"
      },
      "source": [
        "x_data=np.linspace(0,final_pred_y.shape[0], num=final_pred_y.shape[0])\r\n",
        "error = np.array([1.96*(np.std(allres[:,k])/(np.sqrt(allres[:,k].shape[0]))) for k in range(allres.shape[1])]) \r\n",
        "per_05 = final_pred_y + error\r\n",
        "per_95 = final_pred_y - error\r\n",
        "fig = plt.figure()\r\n",
        "params = {'legend.fontsize': 10}\r\n",
        "plt.rcParams.update(params)\r\n",
        "ax = fig.add_subplot(111) \r\n",
        "act = ax.plot(x_data, actual,color = 'red', label='actual')\r\n",
        "pr = ax.plot(x_data, final_pred_y,color = 'blue', label='predicted')\r\n",
        "s1 = ax.plot(x_data, per_05, color='green', alpha=0.4,label = 'Error')\r\n",
        "s2 = ax.plot(x_data, per_95, color='green', alpha=0.4)\r\n",
        "ax.fill_between(x_data, per_05, per_95, facecolor='g', alpha=0.4)\r\n",
        "ax.set_ylabel('PM2.5 Concentration')\r\n",
        "ax.set_xlabel('Time (samples)') \r\n",
        "ax.set_title('Actual vs Predicted')\r\n",
        "ax.legend()\r\n",
        "plt.show()\r\n",
        "plt.close()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YQYp8uNo0q14"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}